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Phase Transition in Strongly Degenerate Hydrogen Plasma 
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Direct fermionic path-integral Monte-Carlo simulations of strongly coupled hydrogen are presented. Our 
results show evidence for the hypothetical plasma phase transition. Its most remarkable manifestation is the 
appearance of metallic droplets which are predicted to be crucial for the electrical conductivity allowing to 
explain the rapid increase observed in recent shock compression measurments. 
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Hydrogen at high pressures remains the subject of 
many investigations, see e.g. [|l], || for an overview. 
Many interesting phenomena, such as a metal-insulator 
transition (MIT) , Mott effect and a plasma phase tran- 
sition (PPT) have been predicted. They occur in situa- 
tions where both, quantum and Coulomb effects are im- 
portant, making a theoretical analysis difficult. Among 
|the most promising theoretical approaches to such sys- 
tems is the path-integral quantum Monte Carlo (PIMC) 
method & |^ which has seen remarkable recent progress, 
e.g. [§[]5)- However, for Fermi systems, these sim- 
ulations are substantially hampered by the so-called 
fermion sign problem. Additional assumptions such as 
fixed node and restricted path concepts have been intro- 
duced to overcome this difficulty §. It can be shown, 
however, that such assumptions do not reproduce the 
correct ideal Fermi gas limit 

Recently, we have presented a new path integral 
representation which avoids additional approximations 
(direct path integral Monte-Carlo, DPIMC) which has 
successfully been applied to strongly coupled hydro- 
gen 10, H, H, see below. In this work we apply the 
DPIMC method to the analysis of dense liquid hydro- 
gen in the region of the hypothetic plasma phase tran- 
sition [|o[ |l|, |ll[ |l^. Computing the equation of state 
land the internal energy, we find clear indications for the 
existence of the PPT - its first confirmation by a first- 
principle method. It is shown that the PPT manifests 
itself by the formation of large metallic droplets which 
are crucial for the plasma transport properties. 

It is well known that the thermodynamic properties 
of a quantum system are fully determined by the parti- 
tion function Z. For a binary mixture of Ng electrons 
and Ni protons, Z can be written as 



ZiNe,N,,V,p) 



Q{N„N„f3)/Njm 
J2jdqdrp{q,r,a;f3). (1) 



Here, q = {qj,q2,... , Qat, } are the coordinates of 
the protons and and a — {ai , . . . , <Jn^ } and r = 
{ri,... ,r7v^} are the electron spins and coordinates, 
respectively, and /3 = l/ksT. The density matrix p 
in Eq. ( p]) is represented in the common way by a path 
integral [O : 
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where A/3 = /3/(n + I) and = 27rfPAl3/me. Fur- 
ther, r"+^ = r" , a' — a; i. e., the electrons are repre- 
sented by fermionic loops with the coordinates (beads) 
[r] = [r, r*^^^ , . . . , r^"-* , r]. The electron spin gives rise to 
the spin part of the density matrix S, whereas exchange 
effects are accounted for by the permutation operator P 
and the sum over the permutations with parity up. 

Following [|, we use a modified representation (||) 
of the high-temperature density matrices on the r.h.s. 
of Eq. (H) which is suitable for cfhcient direct fermionic 
PIMC simulation of plasmas. With the error of the 
order e ^ {(3Ry)'^x/ + 1) vanishing with a growing 
number of beads, we obtain the approximation 
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where x is the degeneracy parameter and t/*, U[ and 
Uf^ denote the sum of the binary Kelbg potentials 
Ip^ , [isll between protons, electrons at vertex 'T', and 
electrons (vertex "/") and protons, respectively. 



In Eq. (|), 



exp[— 7r|^p'' arises from the ki- 



netic energy part of the density matrix of the electron 
with index p, and we introduced dimensionless distances 
between neighboring vertices on the loop, . . 
Finally, the exchange matrix is given by 



exp 



with 
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where the subscript s denotes the number of electrons 
having the same spin projection. From the above ex- 
pressions (0)~(||) one readily computes the internal en- 
ergy and the equation of state: 



f]E = -{Ne + N,) 



f3 



dlnQ 
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In our simulations we used N,, — Ni = 50 and 
n = 20. To test the MC procedure, we have consid- 
ered a mixture of ideal degenerate electrons and classi- 
cal protons for which the thermodynamic quantities are 
known analytically. The agreement, up to the degen- 
eracy parameter x large as 10, has been very good 
and improved with increasing number of particles 
Further, the method was successfully tested in applica- 
tions to electrons in a harmonic trap ||l^ . For the case 
of interacting electrons and protons in dense hydrogen 
we have previously performed a series of calculations 
over a wide range of the classical coupling parameter 
r and degeneracy x for temperatures T > 10, 000 K. 
The analysis of the results has clearly shown a number 
of interesting phenomena, such as formation and decay 
of hydrogen bound states [|[ |[ ^ , including hydrogen 
atoms, molecules, molecular ions, clusters and futher, 
at high densities, pairing of electrons and ordering of 
protons into a Wigner crystal 1^ . 

In this work we present new results which concen- 
trate on the hypothetical plasma phase transition |jl^ . 
For this purpose, we analyze the plasma properties and 
compute the equation of state and internal energy 
(Q) of dense hydrogen along two isotherms, T = 10"* K 
and 5 • 10^ K. Fig. 1 shows pressure and energy vs. 
density at T = 5 ■ 10^ K. For comparison, we also 
include the results for an ideal plasma. As expected. 



due to Coulomb interaction and bound state formation, 
the nonideal plasma results are below the ideal ones. 
We mention that our results are in a good agreement 
with restricted path integral calculations of Militzer and 
Ceperley (Fig. la contains available data points for a 
slightly higher temperature of 6.25 • 10'' K [^. For 
higher temperatures the agreement is very good |l5| . 
Most importantly, at this temperature the pressure in- 
creases monotonically with density, and at high densities 
a continuous increase of the degree of ionization (Mott 
effect) is found. 

However, at T = 10* K the properties of the hydro- 
gen plasma change qualitatively, cf. Figs. 2,3. While the 
overall trend of the pressure. Fig. 2. a, is still a mono- 
tonic increase, in the density region of 0.1 ... 1.5 g/cc 
the plasma exhibits unusual behavior. Inside this re- 
gion the Monte Carlo simulations do not converge to an 
equilibrium state, the pressure strongly fluctuated and 
reached even negative values. Such a behavior is typi- 
cal for Monte-Carlo simulations of metastable systems. 
Note that no such peculiarities appear for densities be- 
low and above this interval as well as for the isotherm 
T — 5 ■ 10* K and for higher temperatures. 

These facts suggest that our simulations have en- 
countered the plasma phase transition predicted by 
many chemical models of partially ionized hydrogen, e.g. 
|lO| , 0, |l^, |l^. According to these models, this is a first- 
order transition with two coexisting phases of different 
degree of ionization. While canonical Monte Carlo sim- 
ulations do not yield the coexisting phases and the coex- 
istence pressure directly, they allow to analyze in detail 
the actual microscopic particle configurations. A typ- 
ical particle arrangement inside the instability region, 
T = 10* K and p = 0.3346 g/cc, is shown in Fig. 3. Ob- 
viously, the protons arrange themselves into large clus- 
ters (droplets) with the electrons (the piecewise linear 
lines show their closed fermionic path) being fairly de- 
localized over the cluster. This is a clear precurser of 
the metal-like state which is found in the simulations 
for densities above the instability region. 

As mentioned above, the plasma phase transition 
appears in many chemical models in the same density- 
temperature range. However, these simple approaches 
become questionable in the region of pressure ioniza- 
tion and dissociation where a consistent treatment of 
all possible pair interactions, including charge-charge, 
neutral-neutral and charge-neutral, is crucial. Further- 
more, these approaches neglect larger bound aggregates 
such as clusters which our simulations reveal to be cru- 
cial in the metastable region. We mention that indirect 
indications for a phase transition have been found in 
recent density functional studies pOl. In this work the 
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thermodynamic properties of hydrogen in the metaUic 
phase were computed (see data points in Fig. 2) and 
enhanced long-wavelength ion density fluctuations were 
observed as the density was reduced to p = 0.799 g/cc 
(the lowest density explored). This led to unusual be- 
havior of the ion-ion structure factor and the effective 
potential which the authors of Ref. |^ interpreted as 
a possible precursor to an incipient metal-to-insulator 
transition. 

Our simulations suggest that the existence of the 
PPT should have a noticable influence on the trans- 
port properties. In fact, when the density changes from 
0.1 . . . 1.5 g/cc, hydrogen transforms from a neutral into 
a metallic fluid. Accordingly, the electrical conductiv- 
ity should increase rapidly. Indeed, shock compression 
experiments have revealed a dramatic increase of the 
electrical conductivity by 4-5 orders of magnitude in a 
very narrow density range of 0.3-0.5 g/cc, jl^, |l^ . The- 
oretical models so far cannot reproduce this behavior 
correctly, predicting either a too early (hopping conduc- 
tivity in the molecular fluid) or too late (free electron 
conductivity) increase pl| . But as the experimental 
data (black circles and crosses in Fig. 2. a) are located 
right inside the PPT region, one has to take into account 
a third conductivity mechanism ~ charge transport via 
electron hopping between individual metal-lik droplets. 
Obviously, this mechanism will be effective inbetween 
the regions where the two other effects dominate and 
thus should allow for a much better agreement with the 
experiments. 

Finally, we mention that our simulations predict a 
PPT for pure hydrogen plasma only. In contrast, no 
PPT was found for a binary mixture of 25% of helium 
and 75% of hydrogen (cf. Fig. 2). 

In summary, we presented direct path integral Monte 
Carlo simulations of dense fluid hydrogen in the region 
of the metal-insulator transition (MIT). Our results 
show evidence for the plasma phase transition which, to 
the best of our knowledge, is its first prediction by a first 
principle theory. Most importantly, we found clear evi- 
dence for the formation of large metallic droplets which 
are predicted to play a crucial role in transport and op- 
tics in the region of the MIT at low temperatures. In 
further investigations we will focus on a more precise 
analysis of the MIT and the plasma phase transition, 
including determination of its critical point and of the 
transport and optical properties of the droplets. 

We are greatful to W. Ebeling, D. Kremp, and 
W. D. Kraeft and S. Trigger for stimulating discussions. 
We also wish to thank D. Ceperley and B. Militzer for 
useful critical remarks. 
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Fig. 1. (a) Pressure and (b) internal energy for 
hydrogen plasma at T = 5 • 10* K vs. density. 1 — 
direct PIMC simulation of this work, 2 — ideal plasma, 
3 ~ restricted PIMC computations at T ^ 6.25 • 10* K 
0- 
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Fig. 2. (a) Pressure (1-5) and electrical conduc- 
tivity (6,7) and (b) internal energy for hydrogen at 
T = 10"^ K vs. density. 1 — direct PIMC simulation 
of this work, 2 — ideal plasma, 3 — direct PIMC 
simulation of a mixture consisting of 25% helium and 
75% hydrogen, 4 — density functional theory [Q, 5 — 
restricted PIMC computations jl^, 6, 7 — electrical 
conductivity of hydrogen (right axis) , 6 — , 7 — . 

Fig. 3. Snapshot of a Monte-Carlo ceU at T = 10"* K 
and p = 0.3346 g/cc. Black circles are protons, dark 
and light broken lines are representations of electrons 
as fermionic loops with different spin projections. 



